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1 INTRODUCTION 

The equation of state (EOS) is one of the most important 
pieces of physics required as input for a stellar evolution 

^1 code. Over a substantial region of the p, T plane the EOS 
^ is relatively simple, and can be computed with an econom- 

j^-^ical code (Eggleton, Faulkner & Flannery 1973; hereinafter 

^^EFF). But along the boundary of this region, at fairly low 
temperature and high density, there are three important pro- 

*TH cesses which were either ignored or treated very crudely in 
|EFF: dissociation of molecular hydrogen, Coulomb interac- 
Q tions, and pressure ionization. A simple contribution for the 
H2 molecule has in fact been included in the code for several 
years, though not published. Coulomb interactions between 
^ the charged particles provide the major non-ideal correction 
to the pressure at the densities and temperatures encoun- 
tered in stars of around a solar mass or less, while they also 
crucially influence the properties of matter at high densities 
and low temperatures. 

In this paper (Section 2) we give a simple but adequate 
approximation to the Coulomb interactions, and we describe 
a somewhat improved, though still very approximate, treat- 
ment of pressure ionization. The result is an EOS which, 
like that of EFF, is completely explicit in its input variables 
(temperature and electron degeneracy), and which is ther- 
modynamically consistent, and very easy to compute. It may 
be reasonably accurate for lower-mass stars (~ 0.1 — 1 Mq), 
as well as for the initial phases of cooling white dwarfs, both 
of which were not very accurately treated by the formulation 
of EFF. It appears to agree reasonably well with the much 
more detailed but time-consuming calculations of Rogers 
& Iglesias (1992), Iglesias, Rogers & Wilson (1992) and 
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We present a simple and efficient, yet reasonably accurate, equation of state, which at 
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Rogers (1994) (hereinafter collectively OPAL), and of Mi- 
halas, Dappen & Hummer (1988) and Dappen et al. (1988) 
(hereinafter collectively MHD). 

We have incorporated the improved EOS into the evolu- 
tion program of Eggleton (1971, 1972, 1973). This program 
has been modified many times in over 20 years, though the 
details of the modifications have not normally been pub- 
lished except occasionally and briefly in relation to specific 
stellar evolutionary problems. Considerable care was taken 
during these modifications to keep the code simple: in fact, 
most changes consisted of removing bits of code that were 
necessary when, for example, double precision was an ex- 
pensive luxury, or 200 Kbytes was the maximum space avail- 
able. A recent paper (Han, Podsiadlowski & Eggleton 1994, 
hereinafter HPE) summarizes the main developments in the 
code, and so we will not repeat them here. However, since 
HPE we have improved not only the EOS, but also the opac- 
ity tables, the neutrino loss rates and the nuclear reaction 
network, and so we describe these changes here, in Section 3. 

In Section 4 we apply the resulting evolution code to 
construct a zero-age main sequence (ZAMS), and evolution- 
ary sequences for a wide range of masses. We discuss these 
results only briefly; we hope to consider various aspects of 
them in a later paper. 



2 THE EQUATION OF STATE 

The main reason for the simplicity of the code of EFF was 
that a Helmholtz free energy F was used in which only 
the free-electron contribution was at all complicated. All 
atomic species, ionized or neutral, contributed only sim- 
ple Maxwell-Boltzmann terms. The electrons alone were al- 
lowed a certain degree of complexity, both as Fermi-Dirac 
particles and as the agents of pressure ionization. This meant 
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first that it was very convenient to use as independent vari- 
ables the temperature T, the electron degeneracy xp (the 
Gibbs free energy divided by kT), and the abundances 
(neutral plus ionized, and so independent of the previous 
two variables); and secondly that the ionization equilibrium 
could be calculated explicitly, since the equilibrium of, say, 
H + + e f> H is not dependent (as, physically, it ought to 
be) on the concentration of, say, He, He + or He ++ . Thus it 
is not necessary to solve simultaneously a complicated and 
highly non-linear set of equations for the ionization equi- 
librium of all species. The inclusion of H2, still as Maxwell- 
Boltzmann particles, only complicates this slightly, requiring 
one to solve a quadratic rather than a linear equation to find 
the overall fraction of electrons that are free. 

The quantities that a stellar evolution algorithm re- 
quires to be calculated from the EOS are only pressure 
p, density p, adiabatic gradient V a , and specific heat at 
constant pressure C p , and also (in some algorithms) their 
derivatives with respect to the independent variables. For 
some purposes, such as the study of oscillations, one may 
also require the compressibility 7 and its derivatives. Of 
course, first derivatives of V a , C p and 7 are in effect third 
derivatives of F, and so we see it as important for thermo- 
dynamic consistency that we start with an expression for 
F which is continuous and sufficiently simple analytically 
that even third derivatives can be extracted without great 
labour. 

2.1 Algorithm 

The formulation of our EOS, like that of MHD, is based 
on the principle of Helmholtz free energy minimization (cf. 
Graboske, Harwood & Rogers 1969, hereinafter GHR). We 
start from a free energy of the form 

F(V,T,N e ,Ni) = -\«T A V 

+ kT ^ Ni {^ Ui V^kT)^ - 1 + w) 

+ F e (N e ,V,T) + AF e (N e ,V,T). (1) 

Here, V is the volume per unit mass (the reciprocal of the 
density p). Other variables and constants have their usual 
meanings, unless specified otherwise later on in this section. 
The first term on the right is due to radiation, and is rela- 
tively trivial to include. The second term is the contribution 
from atomic and molecular particles. The N, are the num- 
bers of particles per unit mass of species i. The only non- 
trivial species that we include are H + , H, H2, and He ++ , 
He + , He. We also include seven heavier elements, C, N, O, 
Ne, Mg, Si and Fe, which are assumed to be fully ionized, at 
all temperatures and densities. Each of these elements con- 
tributes one free electron for every two baryons, except of 
course for Fe which contributes fractionally more. It would 
not in fact be conceptually more difficult (although certainly 
more laborious) to include the ionization of more nuclear 
species. But the effect on the EOS of the is 2 per cent of 
material that is 'metals' is slight. Of course, the effect of the 
state of ionization of the metals on the opacity is enormous, 
but we take opacities in tabular form from work that treats 
the EOS much more rigorously (Section 3). 

The energies X, of ionization or dissociation for 



the six non-trivial species are respectively (in eV) 
0, -13.60, -31.68, 0, -54.40, -78.98 (the zero-point of 
energy taken, for convenience, to be the fully ionized state). 
The partition functions to, are in general functions of tem- 
perature, but except for H2 we approximate these by the 
statistical weights of the ground states: to, = 2 for H and 
He + and Lo t = 1 for the other species. For H2 it is nec- 
essary to take into account the low-energy rotational and 
vibrational modes; we therefore adopt a partition function 
which derives from Vardya's (1960) expression for the pres- 
sure equilibrium constant ifp (H2): 

ch 2 = 6608.8 C(T) ( #^V /2 : e D ^ T -( D ^ T > 2 + ( D ^ T > 3 ,(2) 

where -Dh 2 = 2Xh — Xh 2 = 4.48 eV, the binding energy 
of molecular hydrogen relative to atomic hydrogen, and the 
values of the constants D\, D2, D3 are (in eV) 0.448, 0.1562, 
0.0851. The factor 

C(T) = 1- (l+£^)e- D "> /kT , (3) 

taken from Webbink (1975), truncates the partition function 
which would otherwise diverge for large T. 

The third term of equation (1), F e , is the term that 
treats electrons as pure Fermi-Dirac (FD) particles: 

F e = N e ,pkT-Vp e . (4) 

N e is the number of free electrons per unit mass. Both n e 
(= N e /V) and p e , the electron pressure, are the normal FD 
integrals (see Chandrasekhar 1939), which are explicit func- 
tions of ip, T. We approximate these in the same way as EFF. 
For analytical purposes, the right-hand side of equation (4) 
is to be seen as a function of the independent variables 
N e ,V,T. For computational purposes, however, it is much 
more convenient to take xp,T as independent variables. N e 
has then to be computed from ionization equilibrium (see 
below), so that subsequently V, or equivalently p, is deter- 
mined from the FD integral for n e . EFF gave simple but 
accurate approximations for the three FD integrals, n e ,p e , 
and a third one for the internal energy of the free electrons. 
However, since EFF we have found that greater numerical 
accuracy can be achieved if we evaluate a different integral, 
related to s e = S e /V, where S e is the entropy per unit mass 
of the free electrons (see Appendix A). It should be noted 
that the identities 

are satisfied exactly by these approximations. 

The fourth term of equation (1), AF e , is the term into 
which we endeavour to put all the non-ideal modifications 
which will make the EOS accurate enough for stellar interior 
purposes. Because it involves only N e and not the N,, it is 
guaranteed to be simple to use: the degrees of ionization 
and of dissociation will be explicit functions of xp and T. 
The question is whether we can find a functional form that 
will give sufficient accuracy. We address this in the next 
subsection. 

From F are obtained the pressure, entropy and chemical 
potentials, according to the usual rules: 

p = -(—) 

V dV } T,N e ,Ni 
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2.2 Non-ideal corrections 

The fourth term in equation (1), AF e , contains the contri- 
bution to the EOS owing to the non-ideal effects of Coulomb 
interactions and pressure ionization. As was demonstrated 
above, the simplicity and explicit nature of the algorithm 
rely on the chemical potential deriving from this term be- 
ing dependent only on xp (through the electron density n e ) 
and T. We therefore seek a function g(n e ,T) in order to 
approximate both these effects with a AF e of the form 



AF e = -N e kTg(n e ,T), 
which, using equation (9), implies that 
dg 
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To obtain the electron terms that derive from F e , we have 
used the identities (5). The terms d In to, /d In T in equa- 
tion (7) are zero except for H2. 

The degrees of ionization and of molecular dissociation 
are obtained by minimizing the free energy subject to the 
conservation constraints 



2iV H2 + iV H + iV H + = X/ ibi, 
iV He + iV He + + N He++ = y/4m H , 
and 

N e = N H+ + N He+ + 2N He++ + Z/2m u , 



(10) 
(11) 

(12) 



where X, Y, Z (summing to unity) are the usual abundance 
fractions by mass of hydrogen, helium and 'metals'. These 
lead to the conditions, or 'Saha equations': 



2jUh = flH 2 



fJH = He + /U H +: 
HSe = He + fc+i 

and 

^He+ = He+Hlle-H-, 
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(14) 
(15) 
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In equations (13)-(16) and subsequently, we give the AX, 
in eV; thus k in some places should be in eVK -1 rather 
than ergK -1 . In equation (13), we have divided both sides 
by a factor N e , in order that the left-hand sides of all the 
equations (13)-(16) should be explicit functions of tp, T only, 
given that A/j e also has such a form (see below). Thus equa- 
tions (10)-(16) are seven simultaneous equations for the 
seven unknown abundances. They have a particularly sim- 
ple form: eliminating six of them in favour of, say, Nh_, we 
have only a quadratic equation in Nn to solve. All thermo- 
dynamic properties can then be derived using equations (6) 
and (7), and the internal energy follows from U = F + TS. 



2.2.1 Coulomb interactions 

Formulations for Coulomb interactions have been developed 
in the limits of weak and strong interaction. In the weak limit 
(cf. MHD; GHR), one usually defines the plasma interaction 
parameter 
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r 3 kT 



c 2 



J2 t Njzf + iVefle 

T,i N i 



(19) 



where r 3 is a generalized screening length: r 3 = 
(4ttCV J2, N t /kTV)~ 1/2 . Here, and in the rest of this sec- 
tion, the summation i extends over all ions, with charges 
z,. The factor 9 e = dlnrie/dtp corrects for electron degen- 
eracy (it equals 1 in the limit xp — > — 00, and in the limit 
xp — > 00). In this limit, the free energy contribution is 
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■^ j N l kT 



A 



(20) 



where the factor r < 1, defined by GHR, takes into account 
the finite ion size. This approximation is valid for A ^ 0.3, 
and reduces to the Debye-Hiickel approximation (cf. Cox & 
Giuli 1968) for A < 1 and xp < 0. 

In the strong interaction limit, for xp ^> 0, it is cus- 
tomary to define a different interaction parameter for a one- 
component plasma, 



r.kT 



(21) 



where r 8 is the ion- sphere radius r, = (4nN,/3V)- 1/3 . The 
free energy is then the sum of the contributions of the dif- 
ferent plasma components: 



AF C = ^ AFc ' = - X] N ' kT 9( T ' 



(22) 



In the limit T, — > 00, corresponding to a zero-temperature 
lattice structure, g(Ti) ~ 0.91^ (Cox & Giuli 1968). Slattery, 
Doolen & DeWitt (1980; hereinafter SDD) have performed 
Monte Carlo calculations of the excess free energy of a one- 
component plasma in the liquid phase, for 1 < T, < 160, 
and in the solid phase, for 160 < T, < 300, and given an 
analytical expression for g(Ti) as a fit to their results. 

Although the free energy term AFc clearly depends 
on all ion numbers N,, we can arrive at an expres- 
sion of the form (17) by making two approximations. 
First, we replace the ratios T,,N t /N e and T,,N,zf /N e by 
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No/N e o and Hj(XjZj /mj)/N e o respectively, where iV e o = 
YljXjZj/mj = 0.5(1 + X)/toh is the total number of elec- 
trons (bound or free) and No = YljXj/mj is the total num- 
ber of nuclei (atomic or ionic). The summation here extends 
over all nuclear species j (with charge Zj and mass mj), so 
that these quantities are independent of xp and T. Although 
quantities like Y> t N t may vary greatly, their ratios will do so 
much less, so that this approximation should be reasonably 
accurate even in regions of (partial) de-ionization. Secondly, 
we write the sum Y> t N t g(T t ) as Y> t N t ■ g(T), where T is an 
average plasma interaction parameter. The usual way of av- 
eraging would be r = l£,N,r,/l£,N,, but we find it more 
convenient to use T = C 2 e 2 /r kT = (r 3 /r )A = (A/V3) 2/3 , 
where ro = (47rE 8 Aj/31^) -1 / 3 is the average inter-ion dis- 
tance. 

With these approximations, the explicit expression for 

r is 

T= {-) -J^{j^) { iVeO + e ')>W 

and we model Coulomb interactions with a term in AF e of 
the form 



AF c = -N e kT^g c (T), 



ffc(r) 



[(r/r + a 3 y/^ + ( ai W3/f y/^ 



(24) 
(25) 



This term has the correct asymptotic form in the limits of 
both large and small values of T. The values of the constants 
ai, a,2 and 03 in ffc(r) are found by fitting this expression 
to equation (20) for 0.01 < V < 0.3 and to the expression 
for g(T) given by SDD for 1 < T < 160. The values ai = 
0.897 52 (from SDD), a 2 = 0.768 and a 3 = 0.208 give a 
correspondence better than 4 per cent over the whole region. 
The specific heat C p , which involves the second derivative 
of g c and is the quantity most sensitive to the Coulomb 
term, corresponds within 2 per cent to the values derived 
from SDD. We are therefore confident that this treatment 
is fairly accurate up to T = 160, which corresponds to the 
phase transition to the crystalline state. 



2.2.2 Pressure ionization 

In order to model pressure ionization, we need to choose a 
AF e , and consequential A/j e (equation 9), which will op- 
pose the tendency in equation (14) for hydrogen to become 
steadily more de-ionized as xp increases at constant T. The 
behaviour we want is that hydrogen should start to ionize 
again once the density becomes so large that atoms are sep- 
arated by is do, the Bohr radius. We therefore seek a simple 
term with the properties that, 

(a) for n e <C g>o~ 3 j the term is negligible; 

(b) for n e ~ a^~ 3 , Afj e should roughly cancel the standard 
term xp + 13.60/&T; 

(c) for n e ^> a^ 3 , the term should become large and nega- 
tive, but increasing sufficiently slowly with increasing xp that 
not just N e and n e but also p = n e /N e increase monotoni- 
cally. 

The last condition is required to avoid an awkward multi- 
valuedness in the EOS. If the pressure-ionization term were 



allowed to increase the degree of ionization, and hence N e , 
too rapidly, then p, after increasing with increasing xp, would 
decrease temporarily, before increasing again once the hy- 
drogen is fully ionized. Such multivaluedness would imply 
some kind of phase transition, but this does not appear to 
happen with the much more detailed EOSs of MHD and 
OPAL. We cannot claim, however, to have been very suc- 
cessful in fulfilling (c): our best formulation for (a) and (b) 
still leads to multivaluedness, compressed into a small region 
of the p, T plane, as noted in Section 2.3. 

The following term, containing a number of constants 
(ci . . . C4) which can be varied experimentally to give a good 
fit, appears to be adequate: 

AFpi = -NekTg PI (x,y) (x = n e m H , y = 13.60/&T), (26) 

g PI (x, y) = e" (ci/x)C2 [y + xp + c 3 ln(l + x/c 4 )] . (27) 

The exponential at the beginning of equation (27) takes care 
of point (a), for suitable ci,C2. When x is sufficiently large 
that the exponential is effectively unity, the terms y and xp 
take care of point (b), provided that the logarithmic term 
following is still small. Finally, when x is large enough for the 
logarithmic term to become important, point (c) is nearly, 
though not entirely, satisfied for suitable C3 , C4 . The choice 
of values ci, 02,03,04 = 3, 0.25, 2, 0.03 (ci and C4 in units of 
gem -3 ) provides reasonable agreement with the location of 
ionization zones according to MHD. 

The expression (27) was chosen to ensure that the hy- 
drogen becomes fully pressure-ionized at high density. How- 
ever, at somewhat higher temperature and density it also, 
somewhat crudely, takes care of the helium. The formula- 
tion might be improved by adding further terms like (27) 
which use the helium ionization energy (54.40 eV) instead 
of the hydrogen value (13.60 eV), but such elaboration 
appears unnecessary. Furthermore, we do not separately 
treat the pressure dissociation of H2 which should set in 
at p ^ 10 -2 gem -3 . Therefore the ratio Nh 2 /Nh continues 
to increase towards higher density, but in this region of p, T 
the validity of our EOS becomes questionable anyway. How- 
ever, at high density where hydrogen is fully ionized, the H2 
abundance is also negligible. 

When the material is fully ionized at high pressure, the 
equation of state in effect becomes simple again: we expect 
to get the same EOS as we would get by simply making the 
assumption that ionization is total. However, our term (26), 
(27) affects not only the ionization via A/j e , but also the 
pressure and entropy via equations (6) and (7). This contri- 
bution does not go to zero, as it should, at full ionization, 
unless we subtract from the free energy a compensating term 
thus: 



AF e = AFe + AFpi(iVe, V, T) - AF Pl (N e0 , V, T), 



(28) 



where A^o is the total number of electrons, bound or free, as 
defined above. Since the extra term does not depend explic- 
itly on N e , it does not affect A/j e , but it ensures that the 
effect of pressure ionization on both pressure and entropy 
goes to zero as N e N e o- Of course, the Coulomb screening 
term continues to contribute to both pressure and entropy. 
Hence equation (28), along with equations (23)-(27), is the 
extra term which we incorporate in equation (1). 

Our use of xp in equation (27) slightly conflicts with our 
need to subtract the last term in equation (28), because, 




Figure 1. The equation of state in the logp,logT diagram, for an element mixture with X = 0.7 and Z = 0.02. Dashed lines indicate 
where radiation pressure dominates over gas pressure (to the left of the line p g = p r ) . and where electron degeneracy is important (to the 
right of the line tp = 0), Dash-dotted lines show where the plasma-interaction parameter V (equation 23) equals 1 and 160. The dotted 
lines indicate regions where the pressure term due to Coulomb interactions contributes to the ideal pressure (see text) by more than 0.01 
and 0.1 fractionally. The thin solid lines similarly indicate the regions where the pressure term due to pressure ionization contributes 
to the ideal pressure by more than 0.01, 0.1, and 1. Shaded areas show regions where the ratios jVh 2 /-^h , -^H /-/V H + , N-R e /N-^ e + ' an< ^ 
A^ He +/iV He ++ lie in the range 0.1 to 10. Shown as thick solid lines are a few ZAMS models of stellar interiors of different masses, 
calculated with the updated stellar evolution code (Section 4). 
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whereas we have (in fact, as an input variable) the xp that is 
consistent, via the FD integral, with the quantities N e ,V,T 
of the second-last term in equation (28), we can only ob- 
tain the corresponding tpo for the last term if we invert the 
corresponding FD integral for n e o as a function of tpo , T. 
This would be laborious, and not really necessary since equa- 
tion (27) is very ad hoc. Consequently, we use (but only in 
this section of the code) a much weaker approximation to 
the FD integral for both tp(n e ,T) and tpo(n e o,T): 

^ = 2\pj; f = 1.076 54 a;y 3/2 (l + 0.613 15 a;y 3/2 ) 1/3 . (29) 

This expression approximates xp in the region of large de- 
generacy, where it is important in the above context. 



2.3 Selected results and discussion 

Fig. 1 summarizes the effects of the various contributions to 
the EOS. Shown as a function of p, T are the regions where 
radiation pressure, electron degeneracy and Coulomb inter- 
actions influence the EOS. The dotted lines show the regions 
where Coulomb interactions contribute to the ideal pressure 
(i.e. the pressure obtained without the term AF e ) by more 
than 1 and 10 per cent. This contribution is always negative, 
but never leads to a negative pressure. Note that, because 
of degenerate electron pressure, the relative pressure contri- 
bution decreases again for large enough density, although it 
continues to increase in absolute value. However, for T ^ 1 
the thermal properties of the gas become increasingly dom- 
inated by the Coulomb term - in particular C p increases 
by more than a factor of 2 as T goes from 1 to 160. For 
still larger values of T, crystallization sets in, while for tem- 
peratures kT is fiwp, where lo p is the plasma frequency, a 
quantum-mechanical treatment of lattice vibrations is nec- 
essary, and our EOS becomes unreliable. A simple model 
accounting for these effects would in fact not be difficult 
to implement, and is potentially interesting for the study 
of cooling white dwarfs, but this might also require a bet- 
ter treatment of pressure ionization in the envelopes. Note 
that, for a realistic C-O mixture, the line T = 160 lies at 
A log T fa 1.2 higher than for the H-He mixture shown here. 

The shaded areas in Fig. 1 are regions of par- 
tial ionization and dissociation, namely regions where 
iV H2 /iVH,iVH/iVH+,^He/iVHe+, and iV He+ /iV He++ lie in the 
range 0.1 to 10. The thin solid lines show regions where our 
treatment of pressure ionization modifies the ideal pressure 
by more than 1, 10 and 100 per cent (always in a posi- 
tive sense). Although any prescription for pressure ioniza- 
tion necessarily leads to pressure and entropy contributions, 
a comparison with the OPAL EOS (see below) suggests that, 
in our case, they are largely spurious. Pressure ionization oc- 
curs at roughly the same density for all species, whereas, in 
the much more detailed EOS of MHD, He+, for instance, 
dissociates at higher density than H. For log T is 4.5, pres- 
sure ionization sets in so abruptly that a density inversion 
cannot be avoided. This leads to a discontinuity in the p, T 
plane, shown as the very thick solid line in Fig. 1. Around 
this region, and down to densities J; 10 -2 gem -3 , the va- 
lidity of our EOS is very dubious, although it remains ther- 
modynamically consistent everywhere. As a curious aside, 
the failure of our algorithm to pressure-dissociate H2 partly 
compensates for the spurious pressure term in this region! 



In Figs 2 and 3 we compare values that we obtain from 
the present code with values from the OPAL EOS (made 
available through anonymous ftp by F. Rogers). Fig. 2 shows 
contours of the difference between our pressure and that of 
OPAL, as a fraction of our pressure. It can be seen that 
the agreement is better than 2 per cent over most of the 
region where data are available, although our pressure seems 
to be systematically too large. Part of this difference may 
be attributed to the assumption that all 'metals' are fully 
ionized, although this error can be at most of order Z/2, 
where Z denotes metallicity, or 1 per cent for the assumed 
solar mixture. The largest deviation (up to 5 per cent) seems 
to occur in the region where our pressure ionization term 
contributes most strongly to the pressure. Fig. 3 similarly 
shows differences between the V a of our EOS and that of 
OPAL. These values also agree to better than 1 per cent over 
most of the region, although a systematic deviation of order 
0.04 (about 10 per cent) occurs at low temperatures, just 
below the hydrogen ionization band. It appears to be due to 
the fact that Nb_/N- h + increases too rapidly towards lower 
T, or larger p. This may again be (partly) attributed to the 
assumption of fully ionized metals: supposing that Nb_/N- h + 
is correct as a function of n e ,T, then, if N e is larger than 
it ought to be, p = n e /N e is correspondingly too small, and 
hence Nb_/N- h + will be too large at a fixed p. 



3 OPACITIES, NEUTRINO LOSSES, AND 
NUCLEAR REACTION NETWORK 

Entirely separately from the EOS, we incorporate tables 
of radiative opacity which derive from OPAL, and from 
Alexander (1994, private communication; see also Alexan- 
der & Ferguson 1994a, 1994b) for low temperatures where 
molecular opacity is important. Their values overlap for 
log T between 3.8 and 4.1, and agree to within a few per 
cent in this region. We adopt the OPAL values at log T = 4.1 
and the Alexander values at log T = 3.8, and at intermedi- 
ate temperatures we take values weighted linearly in log T. 
At high temperatures and low densities that fall outside the 
boundaries of these tables, we assume simple electron scat- 
tering (corrected for relativistic effects). Degenerate elec- 
tron conductivities are taken from Itoh et al. (1983), supple- 
mented by values from Hubbard & Lampe (1969) in the par- 
tially degenerate region. From the radiative and conductive 
opacities we construct (by reciprocal addition) a rectangular 
table in log p, log T space, the first incrementing in steps of 
0.25 from —12 to 10, and the second in steps of 0.05 from 3.0 
to 9.3. We use quadratic interpolation in the OPAL data to 
achieve this. Outside the subset of this area for which values 
are available, our values are obtained by extrapolation; such 
values of course have no validity, but we believe that a rect- 
angular, equal-interval mesh is much the most convenient 
format for a stellar evolution code. Our mesh is sufficiently 
fine that we can, during an evolution run, use only linear 
interpolation for a sufficiently accurate model. Of course, 
for such purposes as stability analysis, one needs a more so- 
phisticated approach in order to obtain derivatives that are 
reliable. The code does not note whether some point in the 
star has evolved into the invalid extrapolated area, but, since 
the valid area does in practice cover most of the region where 
stellar interiors are expected to lie, this is not a big problem. 
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Figure 2. The fractional difference (p g — p g) oPAL )/Pg between the gas pressure obtained with our EOS and that of OPAL, as a function 
of log p, log T. The grey-scale is proportional to this quantity, and contours are drawn at the indicated absolute values; solid lines are 
positive or zero values and dashed lines are negative values. The composition is as in Fig. 1. The thick solid line is a 1-Mq ZAMS model. 
The crenellated boundary towards the lower right is the limit of values supplied by OPAL. 



It is in the deep interiors of low-mass main-sequence stars 
that one is most likely to pick up values of opacity outside 
the valid area. Somewhat fortuitously, at the lowest masses 
such interiors are entirely convective, and thus the interior 
opacity (but not the surface opacity) is largely irrelevant. 

Opacity values are now available for a considerable 
range of metallicity Z, as well as at rather smaller intervals 
of hydrogen abundance X than heretofore. In any one evolu- 
tionary run it is normally only necessary to have an opacity 
table for a single zero-age metallicity Zo . For each Zo , we use 
a single abundance variable to interpolate in composition; 
this variable is X in regions where X > and Y otherwise. 
We make use of the fact that, during helium burning, there 
exists a fairly unique relationship between the abundances 
of He, C and O, more or less independent of stellar mass. 



For each Zo we have tables for X = 0.70, 0.35, 0.10, 0.03 and 
0, and for Y = 0.5 — Zo, 0.2 — Zo and 0, and we interpolate 
linearly in X or 7. 

Neutrino loss rates are taken from a series of papers by 
N. Itoh and collaborators. Rates for the photo- and pair- 
neutrino processes are taken from Itoh et al. (1989) and for 
the plasma process from Itoh et al. (1992). These rates are 
all independent of the chemical composition. Rates for im- 
pair bremsstrahlung are taken from Itoh & Kohyama (1983) 
for the degenerate liquid region, and from Munakata, Ko- 
hyama & Itoh (1987) for the partially degenerate region. 
The bremsstrahlung rates are composition-dependent, but 
for the partially degenerate region this dependence is a sim- 
ple proportionality to (Z 2 /A). We therefore find it most con- 
venient to construct two tables, one for the rate e p of photo-, 
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Figure 3. The absolute difference V a — V a opal between the adiabatic temperature gradient obtained with our EOS and that of OPAL, 
as a function of logp,logT. Only the gas component is taken into account. Legend as for Fig. 2. 



pair- and plasma processes combined, and one for the rate eb 
of bremsstrahlung neutrinos for a single element, which we 
take to be oxygen. The total neutrino loss rate is then taken 
to be e v = e p + e b (Z 2 /A)/4, since Z 2 /A = 4 for le O. The 
value thus obtained only deviates from the correct rate at 
very high density and relatively low temperature, and when 
the 'mean' composition is very different from oxygen. Such 
densities and temperatures are only encountered in the cores 
of very evolved stars which are predominantly composed of 
oxygen, so the error introduced is usually small. The tables 
for e p and eb are constructed on the same grid as for opaci- 
ties, for 7 < log T < 10 and < log p < 10. 

It is a feature of the Eggleton code (HPE) that the com- 
position equations are solved simultaneously with the struc- 
ture (and also simultaneously with the mesh-spacing). Even 
with modern computing power it would be uneconomical 
to follow a large number of nuclear species in this way. We 



therefore limit ourselves to five major species - 1 H, 4 He, 12 C, 
le O and 20 Ne; but in effect two more species, 14 N and 24 Mg, 
are followed, by allowing for baryon conservation in the dif- 
ferent burning regions. We also include the abundances of 
28 Si and 56 Fe, but in the present code these elements are 
assumed to remain inert; however, they can be incorporated 
in an extended nuclear network. Our selection of five major 
species allows us to follow the evolution of stars to quite a 
late stage, in principle up to oxygen burning. The rates of 
20 nuclear reactions, as listed in Table 1, are stored as tables 
for 6 < logT < 10. 

The reaction rates are taken from Caughlan & Fowler 
(1988), except for the 12 C(a,7) 16 O reaction for which we 
use the rate from Caughlan et al. (1985), on W. Fowler's 
recommendation. Where two or more reactions are written 
in a chain, the later reactions are taken to be in transient 
equilibrium with the first, and their rates are not included. 
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Table 1. Nuclear reaction network. 
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The first five reactions, comprising the pp-chain, are also 
assumed to be in transient equilibrium with each other, so 
that the (small) abundances of 3 He and 7 Be do not need to 
be followed explicitly, although they are in effect also fol- 
lowed. The reaction involving 18 F is of course an invention, 
for convenience; in effect, 14 N is taken to burn to 20 Ne rather 
than to 22 Ne, so that we can continue to use just the five 
nuclear species itemized above. Something similar applies to 
the last step in the (0,0) reaction. Other possible paths for 
the (C,C), (C,0) and (0,0) reactions are neglected. The 
enhancement of the reaction rates due to electron screening 
is taken into account with the prescription of Graboske et 
al. (1973). 



4 RESULTS OF EVOLUTIONARY 
CALCULATIONS 

To test the resulting code, we constructed a ZAMS from 
~ 0.08 to 250 M , and evolved 24 models of mass 0.64 - 
125 Mq, in steps of about 0.1 in the decimal log of mass, 
from the ZAMS to carbon ignition - except that the four 
lowest-mass stars left the asymptotic giant branch (AGB) 
before igniting carbon, and the four highest-mass stars broke 
down at an earlier stage, for reasons that are not yet clear. 
A subset of these models is shown in Fig. 4, along with 
the observational data that Andersen (1991) selected as the 
most reliable from eclipsing, double-lined spectroscopic bi- 
naries. We used a mixing-length ratio of 2.0, and zero-age 
composition X = 0.7, Z = 0.02. For the mixture of elements 
contained in Z we adopted the ratios derived from Solar sys- 
tem meteorites by Anders & Grevesse (1989). At 'zero-age', 
we take the abundances to be uniform except that, in stars 
above ~ 1.5 M , 12 C is assumed to have reached equilibrium 
through the CN cycle. The models were evolved at constant 
mass and with the Schwarzschild criterion to determine con- 
vective boundaries. 

We were particularly concerned to test the following: 

(a) whether reasonable models could be achieved near the 



bottom of the main sequence, where non-ideal effects most 
strongly affect the EOS; 

(b) whether a model of IMq, at age 4.5 Gyr, would be 
located closer to the observed position of the Sun than we 
have found using the previous version of the code; 

(c) whether models of ~ 2 — 3 Mq would spend sufficient 
time in the blue subgiant region, above the 'hook' but not 
yet into thermal-time-scale evolution, to account for the sur- 
prisingly large number of observed systems (Andersen 1991) 
in this region; 

(d) whether our first giant branch and AGB tracks for 
stars of ~ 0.8 — 8 Mq would be significantly different from 
those of HPE (who used low-temperature molecular opac- 
ities from Weiss, Keady & Magee 1990), and whether this 
could in turn affect the relation found by HPE between ini- 
tial (ZAMS) masses and the masses of white-dwarf rem- 
nants. 

We hope to give a more detailed discussion of these 
points in a future paper. But, even in a detailed discussion, 
we can only make subjective judgments. Regarding point 
(a), we certainly were able to construct models down to ~ 
0.075 Mq, where hydrogen burning gives out. The EOS did 
not produce obvious anomalies in such stars, although their 
interiors are very close to the region (Fig. 1) where difficult 
physics is to be encountered, and where our EOS is probably 
least reliable. Comparison with observation is very difficult, 
since the bolometric correction is very uncertain at these 
low temperatures. 

Regarding point (b), our 1-Mq model passed very close 
to the position of the Sun (shown in Fig. 4), at an age of 5.3 
Gyr. It is both cooler and fainter than the Sun at age 4.5 
Gyr by 1 and 7 per cent respectively. With the EFF approx- 
imation, and older opacities, the discrepancies were much 
worse, about 5 and 20 per cent. However, the abundance ra- 
tios for the Sun as determined by Anders & Grevesse (1989) 
imply a value of Z = 0.0188 for an assumed X = 0.7. We 
constructed an opacity table for this metallicity (by linear in- 
terpolation in Z between 0.01 and 0.02, the values for which 
OPAL data are available) and evolved a 1-Mq model with 
this metallicity and opacity table, assuming the same mix- 
ing length ratio as before. This model was still cooler and 
fainter than the Sun at 4.5 Gyr, but only by as little as 0.2 
and 0.7 per cent. This is certainly within the accuracy that 
can be expected from either the code or the opacity tables, 
so that we did not attempt to 'improve' our solar model by 
adjusting the mixing length and the hydrogen abundance. 

Regarding point (c), we note that, although at least 
six stars in Andersen's (1991) sample are, according to our 
models, evolved beyond the 'hook' at the end of core hy- 
drogen burning, the length of time spent in the subgiant 
region, beyond the hook but on the left of the Hertzsprung 
gap, is not negligible. Fig. 4 attempts to show evolution on 
a nuclear time-scale as a solid line, and evolution on a ther- 
mal time-scale as a dotted line (with occasionally a dashed 
line showing an intermediate time-scale). The fractions of 
the main-sequence lifetime spent in the solid portion of the 
post-hook track, on the left of the Hertzsprung gap, are 5, 
2 and 0.5 per cent respectively for 2, 4 and 8Mq. Although 
0.5 per cent may seem very small, it is still more, by a fac- 
tor of about 5, than the time spent on the whole of the 
subsequent dotted track. We believe that selection effects 
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Figure 4. Hertzsprung— Russell diagram of the ZAMS and several evolutionary tracks, calculated with the updated stellar evolution 
code. The masses of the models are indicated at the starts of the tracks, in solar units. The solid portions of the tracks indicate where 
evolution is on a relatively slow nuclear time-scale, the dotted parts show evolution on a thermal time-scale, and the dashed parts show 
an intermediate time-scale. The different symbols indicate the positions of binary components with well-determined masses, radii and 
luminosities; data for eclipsing, double-lined spectroscopic binaries from Andersen (1991) and data for low-mass, visual binaries from 
Popper (1980). The position of the Sun is indicated by a solar symbol (©). The crosses ( + ) show the position on the giant branches of 
the 1-, 2- and 4-Mq models where the binding energy of the envelope becomes positive (see text). 



for bright, eclipsing binaries are quite likely to push up the 
representation of systems containing post-hook subgiants, 
perhaps by the necessary factor of three or four. Thus at 
1.8 — 2.5 Mq, where Andersen's six stars lie, we are not con- 
vinced that there is a real discrepancy in frequency of these 
systems, especially since the statistical base is not large. Of 
course, a much more stringent test of our models would come 
from a comparison of individual binaries with evolutionary 
tracks, computed for the determined masses and (if avail- 



able) metallicity. We hope to attempt such a study in the 
future. 

Regarding point (d), although our effective tempera- 
tures for highly evolved red giants are substantially higher 
than in HPE, mainly because the molecular opacities from 
Alexander & Ferguson (1994a, 1994b) are significantly dif- 
ferent from those of Weiss et al. (1990), the points on the 
evolutionary tracks (shown by crosses) where the binding en- 
ergy of the envelope goes through zero occur at very much 
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the same luminosity and core mass. Thus we do not have to 
alter significantly the initial/final mass relation obtained by 
HPE. This relation (which includes no free parameters) was 
found by HPE to give a very good agreement with the distri- 
bution of masses of planetary-nebula nuclei (Jacoby 1989). 
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Table Al. Coefficients for the Fermi— Dirac integral expansion. 
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APPENDIX A: MODIFICATION TO THE 
EVALUATION OF THE FERMI DIRAC 
INTEGRALS 

Since EFF, a slight variation has been introduced to the 
evaluation of one of the Fermi-Dirac integrals (cf. Webbink 
1975). EFF gave an approximation to these integrals that 
has the correct asymptotic forms in all four limiting regions 
(tp ^> 1, xp <C — 1, T ^> 1 and T ^ 1) and is also a reasonable 
fit over the whole range of xp and T. The approximation is 
given in terms of / and g defined by 
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In our revised equation of state we continue to use the ap- 
proximations to p* and P* given by EFF (corresponding, 
within a numerical constant, to n e and p e respectively in the 
notation of this paper). We use their fourth-order expansion 
(table 5 of EFF) so that the error is about 0.3 per cent at 
worst. However, rather than calculating the internal energy 
U* directly, Q* was introduced to avoid problems with the 
cancellation of large numbers when finding the electron en- 
tropy. It is defined in a similar way to p* and P* by EFF: 
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the coefficients Q m n are given in Table Al. As in EFF the 
coefficients were chosen so that the state functions are ther- 
modynamically consistent; they obey Maxwell's relations ex- 
actly. The entropy per unit volume of the free electrons is 
then given by 
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and the internal energy by 
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